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Abstract 
> ■ 

' In this paper we report on numerical studies of formation of singularities for the semilinear 

wave equations with a focusing power nonlinearity uu — Au = u*" in three space dimensions. 
We show that for generic large initial data that lead to singularities, the spatial pattern of 
blowup can be described in terms of linearized perturbations about the fundamental self- 
similar (homogeneous in space) solution. We consider also non-generic initial data which are 
' fine-tuned to the threshold for blowup and identify critical solutions that separate blowup 

y—' ' from dispersal for some values of the exponent p. 



1 Introduction 



' One of the most interesting features of many nonlinear evolution equations is the spontaneous 

onset of singularities in solutions starting from smooth initial data. Such a phenomenon, usually 
called "blowup" , has been a subject of intensive studies in many fields ranging from fluid dynam- 
ics to general relativity. Given a nonlinear evolution equation, the key question is whether or 
' not the blowup can occur for some initial data. Once the existence of blowup is established for a 

^ I particular equation, many further questions come up, such as: When and where does the blowup 

occur? What is the character of blowup and is it universal? What happens at the threshold of 
blowup? 

In this paper we consider these questions for the simplest nonlinear generalization of the free 
wave equation: the semilinear wave equation with the power nonlinearity 

Mft — Au = u'', u — u{t,x), x&R^. (1) 

where p > 1 is an odd integer. The sign of the nonlinear term corresponds to focusing, that 
is it tends to magnify the amplitude of the wave. If u is small this term is negligible and the 
evolution is essentially linear (actually one has scattering for t oo). However, if u is large 
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the dispersive effect of the hnear wave operator may be overcome by the focusing effect of the 
nonUnearity and a singularity can form. In fact, it is known that if the energy 

E[u] ^ J (^u^ + {\/uf - J^uP+'^^ d?x (2) 

is negative, then a singularity must form in a finite time This theorem says only that the 
solution cannot be continued beyond certain time but it gives no information on how the solution 
looks like as it approaches the blowup time. Probably the best way to learn about the character 
of singularities is to look at explicit singular solutions. For equation Q it is easy to see that 



2b +1) 



-I, T>0, (3) 



is the exact solution which blows up as t ^ T. This solution is obtained by neglecting the 
laplacian in and solving the corresponding ordinary differential equation uu = ■ The 
question is how typical this explicit singular behaviour is. There are several ways to approach 
this problem. On the analytical side there are Fuchsian methods developed by Kichenassamy 
P] which allow to construct open sets of initial data which blow up on a prescribed spacelike 
hypersurface with the leading order asymptotic behaviour given by the solution 3; . On the 
heuristic side there are numerical and perturbative methods which, albeit non-rigorous, allow to 
gain a more detailed information about the character of blowup. In this paper we take the latter 
approach. 

Our main goal was to show that the solution wq determines the leading order asympotics of 
blowup for generic large initial data and the spatial pattern of convergence to this solution can 
be described in terms of the least damped eigenmodes of the linearized perturbations about mq. 
We did this in the spherically symmetric case 
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Utt — Urr Ur = (4) 

r 

for three representative values of p = 3, 5, and 7. These values corresponds to three differ- 
ent classes of criticality of equation (^). To see this, notice that equation has the scaling 
symmetry: if u{t, x) is the solution, so is 

UL{t,x) = L"u{t/L,x/L), a^j^. (5) 

Under this transformation the energy scales as 

E[ul]= L^E[u], /3=^7-Y' (6) 

hence equation is subcritical for p = 3 (/? < 0), critical for p = 5 (/? = 0) and supercritical for 
p > 5 {p > 0). Since the energy Q is not positive definite, this distinction is not very important 
as far as the generic character of blowup is concerned, however, as we shall show below, it is 
relevant for understanding the behaviour of solutions at the threshold for blowup. 

The paper is organized as follows. In section 2 we discuss self-similar solutions of equation 
Q and analyze their stability. In section 3 we present the results of numerical simulations of 
blowup and demonstrate the universality of the blowup profile. Finally, in section 4 we discuss 
the behaviour of solutions at the threshold for blowup. 
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2 Self-similar solutions 



In order to set the stage for the discussion of singularity formation we first discuss self-similar 
solutions of equation I^J. As we shall see below these solutions play an important role in the 
process of blowup. By definition, self-similar are invariant under rescaling ((SJ, hence in the 
spherically symmetric case they have the form 

u{t,r) = {T-t)-^U{p), " = ^' P=Y^t' 

where T is a positive constant, clearly allowed by time translation invariance. Note that each 
self-similar solution, if it is regular for t < T, provides an explicit example of a singularity 
developing at r = in finite time T from nonsingular initial data - for this reason we shall refer 
to T as the blowup time. Substituting the ansatz ((21 into equation Q one gets the ordinary 
differential equation for the similarity profile U{p) 



(1 - p^)U" + (^^ - (2 + 2a)f^ U' - a{a + 1)U + C/^ = 0. 



(8) 



We consider this equation inside the past light cone of the blowup point {t — T,r = 0), that is 
in the interval < p < 1. It is easy to see that for any p equation © has the constant solution 



Uoip) 



2(P+1) 



(9) 



Of course, this solution corresponds exactly to the solution uq of the original equation Q). It 
turns out that besides this trivial solution, for some values of p there exist also nontrivial profiles. 
The existence of such solutions can be proved by the shooting technique which goes as follows. 
One first shows that at the both endpoints of the interval < p < 1 there exist one-parameter 
families of local analytic solutions which behave, respectively, as 

C^(p)-c+^^^^c-icp)p2 for (10) 

and 

Uip)-b+^(^^{p-l)bP-^byp-l) for p->l, (11) 

where b and c are free parameters. Having that, the strategy for finding solutions which are 
regular in the interval < p < 1 is simple: one shoots the solution satisfying the initial condition 
llll|l at p = 1 towards the center and adjusts the parameter b so that C/'(0) — 0. Applying this 
technique Bizon and Maison ^ proved existence for a countable set of parameters &„ {n — 
0, 1, . . .) which determine analytic similarity profiles Un for p = 3 and all odd p > 7. The first 
few similarity profiles Un for p — 3 and p = 7 are shown in figures 1 and 2. 

The behaviour of solutions Un{p) outside the past light cone, that is for p > 1, depends on 
p. One can show (see |4j) that for p = 3 all n > solutions become singular outside the past 
light cone, namely 

U{p) ^ for some po > 1. (12) 

Po~ P 

In contrast, for p — 7,9, . . . all solutions C/„ remain regular outside the past light cone. 
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Figure 1: The first three similarity profiles Un{p) for p — 3. The index n counts the number of zeros 
of Un{p) in the interval < /o < 1. When continued beyond the past light cone, these solutions become 
singular at some p > 1. 
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Figure 2: The first three similarity profiles Un{p) for p = 7. In contrast to p = 3, here all profiles are 
monotone and have no zeros. For p ^ oo they decay as U{p) ~ p^^^^ ■ 

For p = b there are no nontrivial self-similar solutions - this can be showed as follows. Consider 
the function 



3(5 -p) 2 
4(p-l) ~ (p-l)2 



p + l ^ ' 



This function was introduced by Kavian and Weissler [5] in their study of equation Q . They 
showed that Q'{p) = for p = 5, hence in this case Q is the first integral of equation Since 
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Q(0) = 0, it follows that Q(l) = 0, from which one gets b = [/(I) = (3/4)1/''. rpj^jg coincides 
with Uo so by uniqueness we conclude that Uq is the only similarity profile. 

In order to determine the role of self-similar solutions in dynamics it is essential to analyze 
their stability. To this end we define the slow time t = — ln(T — t) and rewrite equation Q in 
terms of the new variable U (r, p) defined by 



w(i,r) =e"^[/(T,p), 

We get 

Urr + (1 + 2a)Ur + 2pUrp = (1 - p^)Upp + 



p-1 



(14) 



(2 + 2a)p]Up~a{a + l)U + UP. (15) 



The advantage of this formulation is that self-similar solutions of equation |0J correspond now 
to r-independent solutions of equation (|15|) while the asymptotics of blowup corresponds to the 
behaviour at r ^ oo. In order to determine the linear stability of solutions Un{p) we seek 
solutions of l(T5|l in the form U{t, p) = Un{p) + e^^^{p)- After hnearization we get the quadratic 
eigenvalue problem 



,1 - + ^- - 2,1 + e - 2.Af + l.Ur' - + l, - d + 2o,A - A^] I . „, (16) 

Let us consider first the stability of the constant solution Uo. In this case equation becomes 



p+i 



+ A p 



2{p+l) _ p + 3 
p — 1 p — 1 



Near p = the admissible solution has the formal power series expansion 



(17) 



(18) 



k=0 



with the coefficients satisfying the recurrence relation 



A2+(4fc+£±f) A + 2fc (2fc+£±f) -2£±i 
" 2(fc+l)(2fc + 3) 



(19) 



Since Uk+i/ak — > 1 as fc — ^ oo, the series p8f) diverges for p > 1. In order to pass smoothly 
through the point p — I (so that the eigenfunctions are well-behaved also outside the past light 
cone) we impose the condition that the series truncates at the fc-th term 



Uk 



p-1 



A + 2fc 2k 



P 



p-1 



-2P±1=0. 
p-1 



This yields two infinite sequences of pairs of real eigenvalues 

2{p+l) 



Afc 



1 - 2fc, Xt 



p-1 



2k. 



0,1, 



(20) 



(21) 
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There is exactly one positive eigenvalue Aq = 1. It corresponds to the gauge mode which is due 
to the freedom of choosing the blow-up time T. All the remaining eigenvalues are negative hence 
for any p the solution C/q is linearly stable. This suggests that it can appear as an attractor in 
generic evolution. 

Since we do not know the solutions C/„ with n > in closed form, their spectrum of linear 
perturbations can be computed only numerically. Our numerical calculations indicate that the 
solution Un has n unstable modes (apart from the spurious unstable mode corresponding to the 
change of blowup time). For this reason the solutions with n > are not expected to appear 
in generic evolution. However, as we shall see below, the solution Ui with one unstable mode 
appears as the codimension-one attractor in the evolution of specially prepared initial data. 

3 Blowup profile 

Having learned about the stability of the solution uo we are now prepared to interpret the results 
of numerical simulations. The main goal of these simulations was to determine the asymptotics 
of blowup. We solved equation Q using the method of lines which was fourth order accurate in 
space and time. We found that for sufficiently large initial data the amplitude u{t^ r) becomes 
unbounded in a finite time T for some r = rs- More precisely, we have 

liin (T - t)°'u{t, rs) = a = 

which confirms the expectation that the solution uq determines the leading order asymptotics 
of blowup. In this section we wish to show that if the blowup point is at the center, i.e. rs = 0, 
then the spatial pattern of the developing singularity can be described in terms of the least 
damped eigenmodes about uq. 

Using the results of linear stability analysis we can represent the asymptotic approach to Uq 
for T oo (i.e. i — > T) by the formula 

U{t, p) = Uo + Y, Cke^'^Up) + E^fee'^'efc(p), (23) 

k=l fe=0 

where ^k{p), £,k{p) ^^'^ the eigenmodes corresponding to the eigenvalues A^, A^, respectively and 
Cfe, Ck are the expansion coefficients. Keeping the first two least damped eigenmodes we obtain 
the following expansions in terms of original variables (using the abbreviation 5 = T — t) 

P = 3 

<5u(r, t)=V2 + c,5{l-'^^+ c,S' (^1 - |! + ^) + 0{S% (24) 

p = 5 

V6u{t, r) = (^^y + c,S (l - 5) + (co + c/^{l - ^)) + 0{6% (25) 

P = 7 

^1/3 ^(i^ r)=(Jj'\ c,6 (l - ) + c,S'^' + 0{S'). (26) 



2(P+1) 



(22) 
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We claim that these formulae describe accurately the convergence to the blowup profile inside 
the past light cone of the blowup point (t = T,r = 0). The numerical evidence for this assertion 
is summarized in figures 3 and 4 in the case p = 3 (throughout this section we use the case p = 3 
for illustration - analogous results hold for p = 5 and 7). 





Figure 4: For the same numerical data as in figure 3 we plot the deviation of the rescaled solution 
6u{t,r/5) from the constant solution Uq = \/2 at time 6 = 6.09 x 10^'^ The solid line shows the least 
damped eigenmode ci(l — p^)S with the same coefficient ci as in figure 3. 
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Figure 5: For p — 3 the rescaled solution at two moments of time close to the blowup time is shown to 
collapse to the analytic curve F{z) = ^2/(1 + bz^) with b = C1/V2. 



As shown in figure 4 the formula (|24l) accurately describes the blowup profile for large t (i.e. t 
close to T) not only inside the light cone but even slightly beyond. However, the expansions H24t - 
I26II are expected to break down when the linearization is no longer valid; that is, if r^/(5^ ~ 1/6. 
In this transition region the leading order approximation for any p reads 

u{t,r)^^(a + cidi2'^), (27) 



(5" V S ^ 

with a and a defined as in Q and di2 equal to the coefficient of the quadratic term of the ^1 
eigenfunction. This indicates the parabolic scaling 

u{t,r) ^ ^F{z), z = -^. (28) 

Substituting this ansatz into equation Q and dropping the laplacian (which becomes negligible 
as ^ — » 0) we get the ordinary differential equation 

z^F" + (4a + S)zF' + 4a(a + 1)F - 4FP = 0, (29) 

which has a one-parameter family of regular solutions 

a 

(1 -fFz^y 

Comparing H27|) with H28|l and (|30|l we get the matching condition 

di2 



n^) = , .^2.0. - (30) 



aa 



■Cl, (31) 



which e.g. for p = 3 gives h — ci/\/2. The numerical confirmation of this prediction is shown in 
figure 5. 
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We remark that the above result follows immediately from the Fuchsian analysis which pre- 
dicts the leading order asymptotics on a spacelike blow-up curve T{r) in the form 

Ht,r) ^ ^ (32) 
(r(r) - t) 

The blowup time is defined as T = inf T(r). Assuming that this infimum is attained at r = 0, 
we have T{r) ~ T -I- br^ for some 5 > 0. Inserting this into (|32|) we get 

'u(t,r) — —- 5—^, (33) 

5" (1 + 5^)"' ^ ' 

which reproduces and (|^ . 




Figure 6: For p = 3 we plot the rescaled solution approaching the blowup time for initial data fine-tuned 
to the borderline of blowup at rg = and rs > 0. The solid line shows the fit to the analytic prediction 
G{z) = V2/{l + dz*). 

As the coefficient b becomes negative, the first blowup occurs at rs > 0. By fine-tuning initial 
data to the transition between the blowup at rs — and the blowup at > we set b — which 
means that the first eigenmode in the expansion H23|l is tuned away. For such codimension-one 
initial data the formula (|27|l should be replaced by 

- ^ (^a + C2rf23y) , (34) 

where the coefficient (^23 is equal to the quartic term of the ^2 eigenfuction. This gives another 
scaling (see figure 6) 

where 



r 



^23 



■ C2. 



(35) 
(36) 
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4 Threshold for blowup 

Since solutions of equation Q disperse for small initial data and blow up for large initial data, 
there arises a natural question what happens in between. In the following the boundary between 
initial data that lead to dispersion and initial data that lead to singularity formation will be 
referred to as the threshold for blowup. The determination of the threshold for blowup and 
the corresponding dynamics is of great interest in physical models which predict formation of 
singularities, for example in general relativity. This issue can be studied numerically as follows. 
Consider a one-parameter family of initial data (f>{p) such that the corresponding solutions exist 
globally if the parameter p is small or blow up if the parameter p is large. Then, along the curve 
(pip) there must be a critical value p* (or an interval [Pmin^Pmax]) which separates these two 
scenarios. Given two values p small and p large, it is straightforward (in principle but not always 
in practice) to find p* by bisection. Repeating this procedure for different interpolating families 
of initial data one obtains a set of critical data which by construction belongs to the threshold 
for blowup. Having that, one can look in more detail at the evolution of critical data. The 
precisely critical data cannot be prepared numerically but in practice it is sufficient to follow the 
evolution of marginally critical data. Typically, one finds that the evolution of such data has a 
universal (that is family independent) transient phase during which the solution approaches a 
kind of an intermediate attractor. 



p>p* singularity 




Figure 7: A schematic phase space picture of dynamics at the threshold for blowup. 

The heuristic explanation of this behaviour is sketched in figure 7. According to this picture the 
threshold for blowup is given by the codimension-one stable manifold Ws{u*) of an intermediate 
attractor u*, called the critical solution. The critical initial data corresponding to intersections 
of Ws{u*) with different interpolating one-parameter families of initial data, converge^ along 
Ws(u*) towards the critical solution. The marginally critical data, by continuity, initially remain 

^It should be stressed that for conservative wave equations, such as Q, the convergence (which is due to 
radiation of energy to infinity) is always meant in the local sense. 
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close to Ws(u*) and approach u* for intermediate times but eventually are repelled from its 
vicinity along the one-dimensional unstable manifold. Within this picture the universality of 
marginally critical dynamics in the intermediate asymptotics follows immediately from the fact 
that the same unstable mode dominates the evolution of all solutions. The nature of the critical 
solution itself depends on a model - typically u* is a static or a self-similar solution with exactly 
one unstable mode. 

To apply the numerical strategy outlined above we solved equation Q for various one- 
parameter families of initial conditions which interpolate between small and large initial data. 
The results described below do not depend on the particular choice of the family - for concrete- 
ness we present them for the initial data of the form 



u(0, r) = Ar^ exp 



ut{Q,r) = 0, (37) 



with adjustable amplitude A and fixed parameters a and R. Since the initial data are time 
symmetric, the initial profile splits into ingoing and outgoing waves travelling with approximately 
unit speed. Except for very large initial amplitudes for which the singularity forms very fast, 
before the separation into ingoing and outgoing wave occurs, the evolution of the outgoing wave 
does not affect the singularity formation so we shall ignore it. The behaviour of the ingoing 
wave depends on the amplitude A. For large amplitudes we observe the formation of singularity 
at some rs > in finite time T. As A decreases, the blowup point rs decreases also and 
reaches^ rs — Q for some value Aq. As we keep decreasing the amplitude below Aq we eventually 
reach a critical value A* below which solutions do not blowup up. The asymptotic pattern of 
blowup described in the previous section applies to solutions with amplitudes A* < ^ < ^o- 
The character of the threshold for blowup at A* depends on p so we discuss three values of p 
separately. 

p = 7 : 

In this case we identify the critical solution as the n = 1 self-similar solution 

u^{t,r)^{T-t)-"^ Ui{p). (38) 

The numerical evidence for the criticality of solution ui is presented in figure 8. 

According to the picture of critical behaviour described above, the marginally critical solu- 
tions have the following form in the intermediate asymptotics 

u{t,r) = (T-i)"3C/i(p) + C(A)(T-t)-^i-*a(p) + damped modes, (39) 

where is the single unstable mode about ui with the eigenvalue Ai — 11.6442. A small 
constant C{A)^ which is the only vestige of initial data, quantifies an admixture of the unstable 
mode - for precisely critical data C{A*) = 0. We show in figure 9 that the departure from the 
critical solution proceeds in agreement with equation H39|) . 

^The behaviour of the function rg{A) depends on p. As A Aq from above, the function rs(A) decreases 
continuously to zero for p = 3 but for p > 3 it jumps from some rg > to rg = 0. 
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Figure 8: For p = 7 the time derivative Ut{t,r) of marginally critical solutions is plotted for two 
moments of time during the transient phase of evolution and compared to the theoretical prediction 
dui/dt — {T — t)~^''^{Ui/3 + pU[). The parameter T is the same for both curves. 



0.0008 
0.0007 
0.0006 
<^ 0.0005 
^ 0.0004 
0.0003 
0.0002 



supercrit. data 
C= 2.5*10"^^ 
subcrit. data 



C = -2.8*10 



■31 



0.0001 I ' ' ' ' ' ' ^ 

3.627 3.6275 3.628 3.6285 3.629 3.6295 3.63 3.6305 

t 

Figure 9: We plot it~^(t,0) for the pair of marginally critical solutions corresponding to initial data 
II37I I with A = A* ± 10^''^. Initially these solutions are indistinguishable but eventually they split and 
depart from the critical solution towards blowup and dispersal, respectively. The theoretical curves, 
corresponding to equation 13911 for r = 0, with two different fitted coefficients C are superimposed. 
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We know from section 2 that in this case there are no nontrivial self-similar solutions. However, 
since for p = 5 the energy is scale invariant, static solutions with finite energy are possible. 
Indeed, it is well known that equation |0J has the finite energy solution 

usW = (l + ir2)-5. (40) 

Rescalings of this solution generate the orbit of static solutions Ug = L~^/'^us{r / L). To deter- 
mine the linear stability of this solution we plug u{t, r) = us{r) + e^^*v{r) into Q and linearize. 
We get the eigenvalue problem in the form of the radial Schodinger equation 

-v" ~-v' + Vv = k^v, V = -- . (41) 



(l + ir2)2 



Notice that the perturbation induced by rescaling 

d 



L=l 



(l + ir2) 



satisfies equation (|41|l for fc^. This is so called zero mode. Since the zero mode has one node, 
it follows by the standard result from Sturm-Liouville theory that the potential V has exactly 
one bound state, kf < 0, which means that there is exactly one growing mode e'''i*wi(r), where 
Ai = ^ —k\. Numerical calculation gives Ai « 1.1. 

Thus, according to our preceding discussion, the solution us is a candidate for a critical 
solution. To verify this, Szpak [B] has investigated the nonlinear evolution of the growing mode. 
For initial data of the form u(0, r) — us(r)~\-ev\[r), Mt(0,r) = eAitii(r), he found that, depending 
on the sign of the amplitude e, the solution either disperses or blows up in finite time. This 
confirmed the expectation that in fact us is the critical solution sitting on the saddle separating 
blowup from dispersal. Applying bisection to the family of initial data (|37|l we have obtained 
the solution us as the intermediate attractor with pretty long lifetime. We refer the reader to 
[HI for more details, in particular the analysis of convergence to us- 

p = 3 : 

In this case we were not able to identify a critical solution because of two reasons. First, 
in contrast to the cases described above, in p = 3 we do not have a good candidate for the 
critical solution. The only potential candidate is the self-similar solution ui with one unstable 
mode, however, as mentioned in section 2, this solution is singular outside the past light cone 
of the blowup point and therefore cannot be a bona fide critical solution. Second, we face the 
following difficulty when trying to determine A* . As we approach the expected value of A* from 
above the wave initially shrinks but at some later time ti it bounces back and expands outside 
with decreasing amplitude. During this period of evolution the amplitude of outer wave front 
decreases faster then the amplitude of the solution at the center, so a flat, slowly decreasing with 
time central region develops. After some time t2 this central part of solution returns and starts 
growing again to form a singularity at r = 0. If we decrease A further, the time of bounce ti 
almost does not change but the return time t2 increases significantly. Therefore, approaching A* 
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we have to evolve the solution longer and longer on larger and larger grids. Since the numerical 
grid is always finite, we cannot tell if an expanding wave which leaves the grid represents a 
genuine dispersion or a singular solution with large return time t2- Figures 10 and 11 illustrate 
this difficulty. 




Figure 10: For p = 3 wc plot the snapshots from the evolution of the wave that has bounced Ijacfc from 
the center. After the bounce the amplitude at the center initially decreases but later the wave returns 
and the amplitude starts growing again. 
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Figure 11: The same data as in figure 10. The first minimum of \ I Umax corresponds to the bounce. 
The second local maximum corresponds to the return. 
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5 Conclusions 



We have studied formation of singularities for the spherically symmetric semilinear wave equation 
with the focusing power nonlinearity for three representative values of the exponent p: p = S 
(subcritical case), p = 5 (critical case), and p = 7 (supercritical case). We showed that in 
all these cases the asymptotic behaviour of blowup can be understood in terms of decaying 
perturbations about the fundamental (homogeneous in space) self-similar solution. We showed 
also that the nature of the critical solution, whose codimension-one stable manifold separates 
blowup form dispersal, depends on p: for p = 7 the critical solution is self-similar while for p = 5 
it is static. For p = 3 we were not able to identify a critical solution determining the character 
of the threshold for blowup in this case remains an open problem. 
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